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Abstract 

We investigate the effect of early chemical freeze-out on radial flow, elliptic flow and HBT radii by 
using a fully three dimensional hydrodynamic model. When we take account of the early chemical 
freeze-out, the space-time evolution of temperature in the hadron phase is considerably different 
from the conventional model in which chemical equilibrium is always assumed. As a result, we find 
that radial and elliptic flows are suppressed and that the lifetime and the spatial size of the fluid 
are reduced. We analyze the pt spectrum, the differential elliptic flow, and the HBT radii at the 
RHIC energy by using hydrodynamics with chemically non-equilibrium equation of state. 
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I. INTRODUCTION 



One of the main topics in the physics of relativistic heavy-ion collisions is to give a detailed 
description of space-time evolution for the hot and dense nuclear matter by using a dynam- 
ical model such as kinetic transport models [|TJ or hydrodynamics 0. Recent experimental 
data at the RHIC energy suggest that we may see the significant effect of jet quenching in 
transverse momentum distribution for neutral pions [[| or azimuthal asymmetry for charged 
hadrons Q]. So hydrodynamic simulations of expanding hot and/or dense matter are in- 
dispensable in quantitatively estimating the effect of the medium on the jet quenching ||. 
One of the authors (T.H.) has already built a fully three dimensional hydrodynamic model 
which describes not only central but also non-central collisions @. In contrary to other 
hydrodynamic models, e.g., (2+1) dimensional models with the Bjorken's scaling solution 
|7], [|, or (3+1) dimensional models with cylindrical symmetry along the collision axis 



T0| , PI, [12| , |T3l , full (3+1) dimensional hydrodynamic models [[14], [15], [T^, [17], [18], [TIJ enables 



us to obtain the rapidity dependence of particle distribution, elliptic flow, and HBT radii in 
non-central collisions. 

Statistical and hydrodynamics-motivated models give us the characteristic temperatures 
in relativistic heavy-ion collisions. These temperatures are very useful in understanding what 
happens in collisions. From fitting the model calculation of particle ratios for hadrons to 
the experimental data, the chemical freeze-out temperature T ch is obtained at each collision 
energy [[H], ^T, E2J |23|, |24], 2pJ. On the other hand, we can obtain the thermal (kinetic) freeze- 



out temperature T th from the slopes of transverse momentum distribution by assuming the 
radial flow profile |22, 26], [27|. The temperatures obtained above are usually different from 
each other at the AGS, SPS, and probably RHIC energies: T ch ~ 160-200 MeV while T th ~ 
100-140 MeV. The chemical freeze-out parameters (T ch , fj,^) at various collision energies seem 
to be aligned on one line in the T-fi^ plane. This line is sometimes called "chemical freeze- 
out line" and can be parameterized by the average energy per particle < E > / < N >~ 1 



GeV Those analyses indicates that the system undergoes first the chemical freeze- 

out where the observed particle ratios are fixed and next the thermal freeze-out where the 
shape of the transverse momentum distribution is fixed ||28|| . Since the time scale of the 
hydrodynamic evolution is comparable with (or smaller than) that of the inelastic collisions 
between hadrons, the number changing processes are likely to be out of equilibrium. This 
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is the reason why there are two sequential freeze-out processes in relativistic heavy-ion 
collisions. To give a more realistic description of the temporal and spatial behavior of the 
hot and dense matter, the above pictures should be included in the model. This makes 
us approach the comprehensive understanding of the deconfmed matter, the Quark Gluon 
Plasma (QGP). 

In this paper, we incorporate the different freeze-out temperatures, T ch and T th , into 
hydrodynamics and discuss how the early chemical freeze-out affects the space-time evolu- 
tion of fluids and the particle spectra [p9fl . In the ordinary hydrodynamic calculations, one 
assumes both chemical and thermal equilibrium and consequently T ch = T th which is to 
be determined from comparison of the slopes of transverse spectrum with the experimental 
data. If the system obeys the above picture of early chemical freeze-out, those ordinary 
hydrodynamic models can hardly reproduce the particle ratios due to the smallness of the 
chemical freeze-out temperature. As a result, the number of resonance particles at the ther- 
mal freeze-out becomes too small. The physics at the low transverse momentum is largely 
affected by resonance decays after thermal freeze-out. For example, the low pt enhancement 
in the transverse momentum spectrum for pions can be explained by the contribution from 
resonance decays |3(], [JTJ. The slope of the p t spectrum for pions directly emitted from 
the freeze-out hypersurface is almost constant. On the other hand, data for pions at low pt 
shows a steeper slope than those in medium p t region (~ 1 GeV/c). So one cannot reproduce 
the low pt enhancement only from direct pions. Pions from p, uj, or A decays show steeper 
Pt spectrum than direct pions. We naturally explain the low pt enhancement seen in the 
experimental data by resonance decays. Another example is the reduction of the second 
Fourier coefficient of the azimuthal distribution for pions at low p t ||32|| . From the exact 
treatment of the decay kinematics, pions from resonance decays can pretend out-of-plane 
elliptic flow even when the hydrodynamic flow shows in-plane elliptic flow. This dilutes 
elliptic flow from direct pions. Therefore the early chemical freeze-out must be included in 
hydrodynamics in order to analyze not only the particle ratios but also the single particle 
spectra and the azimuthal asymmetry. 

This paper is organized as follows: We construct the equation of state for the hadron phase 
with and without chemical equilibrium in Sec. II. We parameterize the initial condition for 
hydrodynamic simulations in the full three dimensional Bjorken coordinate in Sec. [TIT]. By 
using these equations of state and the initial condition, we perform hydrodynamic simula- 
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tions in full three dimensional space at the RHIC energy and compare space-time evolutions 
with each other in Sec. [IV]. We analyze the particle spectra and the azimuthal asymmetry at 
the RHIC energy and discuss the effect of early chemical freeze-out on observables in Sec. |VJ. 
We also analyze the two-pion correlation functions to see the effect on the hydrodynamic 
evolution in Sec. |VJ. Summary and discussions are given in Sec. |V1J| . 

II. EQUATIONS OF STATE 

We assume the following picture of space-time evolution for hot and/or dense matter 
produced in relativistic heavy-ion collisions. First, the huge number of secondary partons 
are produced and both chemical and thermal equilibrium among these partons are achieved 
in the early stage of collisions. The initial dominant longitudinal flow and the large pressure 
gradient perpendicular to the collision axis cause the system to expand and cool down. 
When the temperature of the system reaches the critical temperature T c , the hadronization 
starts to occur. Just after the hadronization finishes, the chemical freeze-out happens at 
T ch (< T c ). Below T ch , the ratios of the number of observed particle are fixed. Even after 
chemical freeze-out, the system keeps thermal equilibrium through elastic scattering. Finally, 
all hadrons are thermally frozen at T th (< T ch ). If we neglect dissipation in the space-time 
evolution of nuclear matter, we can apply the hydrodynamic equations for the perfect fluid, 
d^T^ = and c^n^ = 0, where T^ u — (E + P)u^u u — Pg^ u and = h-qu^ are energy 
momentum tensor and baryon density current, respectively. E, P, riB, and are energy 
density, pressure, baryon density, and four fluid velocity. With the help of thermodynamical 
identities and the baryon density conservation, the first equation is rewritten in d^s^ = 0, 
where s M is the entropy current. These equations mean that a fluid element evolves along 
an adiabatic path {u-q/ s = const.) in the T-/xb plane. We assume this fact is approximately 
valid even below chemical freeze-out line. 

When N stable hadrons in the equation of state (EOS) undergo chemical freeze-out across 
the chemical freeze-out line, we can introduce chemical potentials /ij associated with those 
hadrons. Then we may construct the EOS in the (N + 2) dimensional space (N for ^ and 2 
for T and hb)- This causes a serious problem when we numerically simulate a hydrodynamic 
model with finite resources of memory. Since the chemical potential of each hadron depends 
on T ch , /Xg 1 , T, and fiB during expansion along an adiabatic path, we can restrict ourselves 
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to the EOS in a two dimensional hypersurface \ii = /Xi(T, /xb) embedded in the (N + 2) 
dimensional space. Thus we need not prepare such a large dimensional table of the EOS 
anymore. When we obtain /ij at a point (T, //b) below the chemical freeze-out line, we 
need the information at the chemical freeze-out point (T ch ,/ig h ) somewhere on the chemical 
freeze-out line and have to go back to the point along an adiabatic path. Since the adiabatic 
path itself is obtained from the thermodynamical variables and fii(T, //b), we have to solve 
the problem self-consistently. In this paper, we restrict our discussion to the zero baryonic 
chemical potential where the adiabatic path becomes a trivial one, ub/s = 0. This is a good 



approximation in Au + Au ^/snn = 130 GeV collisions in which p/p ~ 0.6 |33], |34|, |35|, |36 
or the resultant chemical freeze-out parameter /ig h ~ 50 MeV p4| . 



We construct three models EOS to compare the space-time evolution of fluids. These 
models describe the first order phase transition between the QGP phase and the hadron 
phase at T c = 170 MeV. We suppose the QGP phase is composed of massless u, d, s 
quarks and gluons and that it is common to three models. The EOS for the QGP phase is 
P = (E — AB)/3, where B is a bag constant specified later. For the hadron phase, we choose 
three different models EOS as follows. 



A. Chemical Equilibrium 

The first model is an ordinary resonance gas model in which complete chemical equilib- 
rium is always assumed (model CE). This model is employed for the sake of comparison with 
the other models. We include strange and non-strange hadrons up to the mass of A(1232). 
Energy density and pressure are as follows: 

E = S^J±- [ d 3 y y/pt + mj 

P = ^E T (|jl/ (|3 P lo g{ lTeXp - (^p' + mf-^/T }, (2) 

where /Zj = due to complete chemical equilibrium in this model. Here the upper and 
lower signs correspond to bosons and fermions. We neglect the excluded volume correction 
which largely affects the hadronic EOS in the high baryon density region. From the Gibbs's 
equilibrium condition at T c = 170 MeV, we obtain the bag constant B 1 ^ = 247.2 MeV and 
the latent heat AE ~1.7 GeV/fm 3 . 
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B. Chemical Freeze Out 



The second model is the simplest one which describes the picture of chemical freeze-out 
(model CFO). Below T ch , we assume the numbers of all hadrons included in the EOS 
are fixed and that the particle number densities obey d^n? = 0. We introduce a chemical 
potential Hi(T) associated with each species so that iV$ becomes a conserved quantity. From 
the conservation of entropy = 0, the ratio of the particle number density to the entropy 
density below the chemical freeze-out temperature obeys 

nj{T, iij) nj(T ch ,Hi = 0) 

s(T, {m}) ~ s(T dl ,{ t i i } = 0) [ ' 

for all hadrons along the adiabatic path. Here we assume T ch = 170 MeV, which is consistent 
with a recent analysis based on a thermal model at the RHIC energy From Eq. (^), we 
obtain a chemical potential as a function of temperature for each hadron. The Hi(T) ensures 
to keep the ratios of the number of each hadron throughout the space-time evolution of a 
fluid element without explicitely solving d^n^ = 0. All chemical potentials are functions 
of temperature, so the thermodynamical variables depend only on temperature even after 
chemical freeze-out. 



C. Partial Chemical Equilibrium 

The third model represents a more realistic EOS than the second one. The following 
model is first discussed in Ref . [ 37J . The observed particle numbers are always composed of 



the contribution from direct particles and resonance decays, i.e., N n = N n + Ylv^di-^Ni. 
Here d is an effective degree of freedom which is a product of the degeneracy d and the 
branching ratio B. So some elastic processes with large cross sections (e.g., tttt — > p — > tttt, 
ttN — > A — >• ttN, tcK — * K* — ► ttK) can be equilibrated even below T ch as long as the 
equality 

nj{T, iij) nj(T ch ,Hi = 0) 

s(T,{^}) ~ s(T ch ,{ f i l } = 0) [ ' 

is kept instead of Eq. @. We regard ir, K, rj, N, A and S as "stable" particles 1 and that all 
chemical potentials can be represented by chemical potentials associated with these stable 



1 Here we use a term "stable" when the lifetime of a hadron is much longer than that of the fluid (~ 10 
fm/c). 
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FIG. 1: Chemical potentials for pions (dashed lines), p mesons (solid lines), and u mesons (dotted 
lines) for the models CFO and PCE. 

particles, e.g., \i p = 2p, n , p,x* = pL n + Hk, A*a = + fJ>N> an d so on - Thus the third model 
describes the partial chemical equilibrium (model PCE) even below T ch 131. It should 



be noted that, after chemical freeze-out, hn = Hn (¥" 0) with keeping baryonic chemical 
potential \ib = in our model. 

The model PCE employed here is not the only one to describe the partial chemical equilib- 
rium. There may be other choices for stable particles or other processes to be equilibrated 
in this model. Various models should be checked by future precise experimental data of 
particle ratios. 



D. Chemical Potentials and Equations of State 

Figure [I] shows the chemical potentials for n, p, and u mesons for the models CFO and 
PCE. The difference of chemical potentials between hadrons depends only on its mass in 
the model CFO, so p, u (T) behaves like fi p {T) due to the small mass difference. Both are 
almost linearly increasing with decreasing temperature. On the other hand, these chemical 
potentials differ from each other in the model PCE. It results from each elementary process, 
i.e., p <-> mi {pi p = 2pt n ) and uj <-> txtxtx = 3 x 0.88//^), where the branching ratios are 
from Ref . |3T3[ . 

One can easily evaluate EOS's for these models by inserting chemical potentials obtained 
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FIG. 2: (a) Pressure as a function of energy density, (b) Temperature as a function of energy 
density. The dashed, dotted, and solid lines correspond to the models CE, CFO, and PCE. 

above into Eqs. (|]) and @). We represent pressure and temperature as functions of energy 
density for three models in Fig. |[ We find in Fig. § (a) that pressure as a function of energy 
density is similar to each other. On the other hand, temperature as a function of energy 
density in the models CFO or PCE in Fig. [| (b) is deviated from the model CE. Since the 
resonance population of the models CFO or PCE is larger than that of the model CE due 
to the chemical freeze-out, the energy density at a fixed temperature in the hadron phase is 
also large in those models. Conversely, temperature in the models CFO or PCE at a fixed 
energy density is smaller than in the model CE. This fact is very important in qualitatively 
understanding the difference of the space-time evolution of fluids among three models as 
shown in Sec. [TV]. 
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III. INITIAL CONDITIONS 



We numerically solve the hydrodynamic equation in the full three dimensional Bjorken 
coordinate (r, r] 8 , x, and y) which is relevant to analyze heavy-ion collisions at the collider 
energies. Here r = \Jt 2 — z 2 and r] s = (1/2) log[(t + z)/(t — z)\ are the proper time and the 
space-time rapidity, respectively. The x axis is parallel to the impact parameter vector and 
the y axis is perpendicular to the x axis in the transverse plane. As one of the authors (T.H.) 
has already pointed out ||, the main reason to employ the Bjorken coordinate rather than 
the Cartesian coordinate is a practical one: This considerably reduces numerical efforts such 
as the long lifetime of fluids (~ 100 fm/c) and less numerical accuracy near the light cone. In 
addition to this reason, there is also a physical reason why we avoid to employ the Cartesian 
coordinate. When one simulates the space-time evolution in the Cartesian coordinate, one 
gives the initial condition at a constant to- If is regarded as the thermalization time, this 
initial condition implies that the thermalization occurs first from the forward (backward) 
space-time rapidity in T-rj s coordinate. 2 This is somewhat unrealistic because the multiplicity 
in the forward rapidity region is much smaller than the one at midrapidity. This may cause 
a crucial problem when one discusses the rapidity dependence of radial and elliptic flows at 
the collider energies, since they are sensitive to the thermalization of the system. 

We choose initial conditions in the Bjorken coordinate so as to reproduce pseudorapidity 
distribution in Au+Au 130A GeV central (0-6 %) collisions obtained by the PHOBOS 
Collaboration [[HJ. At the initial time To = 0.6 fm/c, the initial energy density for central 
collisions can be factorized as 

E(x,y,r} s ) = E max W(x,y;b)H(ri s ). (5) 

We assume the transverse profile function W(x,y; b) scales with the impact parameter b in 
proportion to the number of binary collisions 



W(x,y;b) oc T + T_, T± = T(x±b/2,y), (6) 
where T(x, y) is a thickness function with the standard Woods-Saxon parameterization for 



2 Supposing to = 1 fm/c, the corresponding initial proper time becomes To = io/7JV ~ 0.01 fm/c near the 
beam rapidity region at RHIC energies. 
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nuclear density 

T(x,y) = 



dzp{x,y,z), p{x,y,z) 



Po 



(7) 



exp[( v /x 2 + y 2 + z 2 - Rq)/5] + 1 ' 

Here p — 0.17 fm -3 is the saturation density, 5 = 0.54 fm is the diffuseness parameter, 
and Ro = 1.12A 1 / 3 — 0.86A -1 / 3 fm is the nuclear radius. The proportional constant in 
Eq. (H) is fixed from the condition W(0, 0; 0) = 1. The longitudinal profile function H(rj B ) is 
characterized by two parts || [13], |42fl : it is flat near r] s ~ and smoothly connects to vacuum 
as a half part of the Gauss function in the forward and backward space-time rapidity regions 

(I Vs -Tfeo I -%at/2) 2 



cxp 



2 ? ?Gauss 



Vs - VsO I -%at/2) 



(8) 



The length of a flat region r/fl at and the width of the Gauss function r/Gauss are adjustable 
parameters to be determined by the experimental data, especially (pseudo)rapidity distribu- 
tion. In symmetric collisions with the vanishing impact parameter, we expect the symmetry 
E(x,y, —rj s ) = E(x,y,r] s ) holds at the initial time. On the other hand, in non-central (or 
asymmetric) collisions, we can shift the energy density by rj s o which is identified with the 
center of rapidity for each transverse coordinate 



r] s0 (x,y;b) = - log 



(TL + T + )7at + (T_ - T + )j N v N 



(9) 



_(T_ + T + )>y N - (T_ - T + )>y N v N \ ' 

where vn and 7jv are, respectively, velocity and Lorentz 7 factor of an incident nucleon in 
the center-of-mass system. For illustrations of the initial energy density, see Ref. 0. The 
initial longitudinal flow is the Bjorken's solution [[HJ], i.e., the fluid rapidity ^(to) is equal 
to the space-time rapidity rj s . It should be noted that this is merely an initial condition 
and that Y { ^ rj s after initial time due to the pressure gradient directed to the 7/ s axis. The 
transverse velocities vanish at r and are to be generated only by the transverse pressure 
gradient. 

Initial parameters in hydrodynamic simulations are so chosen as follows: -E max — 35 
GeV/fm 3 , r/fi a t = 5.8, r/Gauss = 0.2, and b = 2.4 fm. These values lead us to reproduce the 
pseudorapidity distribution in Au + Au central collisions at ^/snn = 130 GeV observed 



by the PHOBOS collaboration [40|. The parameters are adjusted for the model PCE with 
T th = 140 MeV. We also use the same values for the other models EOS for the sake of 
comparison, although this causes the pseudorapidity distribution slightly deviated from the 
experimental data. 
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IV. TIME EVOLUTION OF FLUIDS 



We simulate the space-time evolution of the fluid in the full three dimensional space 
with the initial conditions and the EOS's discussed in the previous sections. First, we pick 
up a fluid element at the central point (x = y = r] s = 0) and pursue its time evolution till its 
temperature reaches T = 100 MeV. Figure [3] shows the time evolution of (a) energy density 
and (b) temperature at the center of the fluid for three models EOS. As far as the time 
evolution of energy density, we cannot distinguish each other. This is easily understood by 
Fig. ^| (a): Energy density evolution is completely governed by the EOS, i.e., P(E) and the 
three models are very similar to each other. We transform these results from energy density 
to temperature by using Fig. |2] (b). For the time evolution of temperature, chemical freeze- 
out makes a substantial difference between the model CE and the model CFO or PCE. If 
we suppose thermal freeze-out occurs at the constant temperature, the system in which the 
property of the chemical freeze-out is considered has a fate to be thermally frozen earlier 
than the conventional model CE. The early chemical freeze-out makes the hadron phase cool 
down more rapidly |45j . 

Next, we show how the early chemical freeze-out affects the spatial size of the fluid. 
Figure ^ represents the time evolution of freeze-out hypersurfaces at y = r] s = for (a) the 
model CE and (b) the model PCE. Here the hypersurfaces in Fig. 4 correspond to various 
T th = 100, 120, 140, and 160 MeV which are within a plausible range for thermal freeze-out 
temperature. We find that the early chemical freeze-out reduces not only the lifetime of 
the fluid but also its spatial size and that the fluid does not expand so explosively for the 
model PCE. Since the two-particle correlation function is sensitive to the spatial size and 
the lifetime of the fluid, it is interesting to see the effect of early chemical freeze-out on the 



HBT radii. Detailed analyses of pion interferometry will be discussed in Sec. [VI . 

Figure [| shows the thermal freeze-out temperature dependences of average radial flow at 
midrapidity < v r >|^ s =o, where v r = ^v"l + v%. Radial flow is generated by the pressure 
gradient, so it contains informations about the EOS. From this figure, the radial flow is 
suppressed when we take account of the early chemical freeze-out. At T th = 140 (120) MeV, 
the average radial flow for the model PCE is reduced by 17.7 (22.5) % from the one for the 
conventional model CE. 

It should be noted that the difference between the thermal freeze-out temperature T th 
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FIG. 3: (a) Time evolution of energy density at the center of the fluid, (b) Time evolution of 
temperature at the center of the fluid. The dashed, dotted, and solid lines correspond to the 
models CE, CFO, and PCE. 

and the thermal freeze-out energy density E th plays an important role in analyses of particle 
spectra. It may be claimed that there are no significant differences among three models when 
one shows the results in Figs. ^ and [5] as functions of the thermal freeze-out energy density 
E th , not temperature T th . It is indeed true when one discusses only the hydrodynamic 
behavior. The shape of particle distribution, especially pt spectrum, are determined by 
the thermal freeze-out temperature T th and the flow in the hydrodynamic model. So this 
difference is clearly meaningful when we compare the numerical results of particle spectra 
with experimental data as shown in the next section. 
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FIG. 4: Space-time evolution of thermal freeze-out hypersurface for (a) the model CE and (b) the 
model PCE. 
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FIG. 5: Thermal freeze-out temperature dependence of average radial flow. The dashed, dotted, 
and solid lines correspond to the models CE, CFO, and PCE. 
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V. SINGLE PARTICLE SPECTRA AND AZIMUTHAL ASYMMETRY 



Hydrodynamics has a limited prediction power, so the calculation is really meaningful 
only after tuning the initial parameters and reproducing the single particle spectra. The mo- 
mentum distribution for particles directly emitted from a (thermal) freeze-out hypersurface 



can be calculated through the Cooper- Frye formula [JE] 

dN_ _ d f ydo 

d*p (27r)3y s exp[(p-u-/i i )/T th ] T l- 1 ' 

Here, £ and da^ 1 are the thermal freeze-out hypersurface and its element. w M is the four 
fluid velocity. — (+) sign is for boson (fermion) and d is the degeneracy of particles under 
consideration. This formula merely counts the net particles passing through the hypersurface 
X rather than decoupling from the system. Although it has a problem on the negative 
number in the treatment of time-like freeze-out hypersurface |[37| , this is widely used in 
almost all hydrodynamic models. The observed spectra always contain the contribution 
from resonance decays. We assume all of the resonance particles in the EOS are also emitted 
from a freeze-out hypersurface and that they decay into stable particles. Taking account of 
the decay kinematics, we easily obtain the single particle spectra from resonance particles 
|3T| , |32|. 3 For further details to calculate the spectra from resonance decays, see Appendix. 
Since the results from the model CFO is similar to the ones from the model PCE, we hereafter 
concentrate our discussions on the models CE and PCE. 



A. Spectra and Flow for Charged Hadrons 

Figure || shows the pseudorapidity distribution of charged particles in Au + Au 130A 
GeV collisions. We choose T th = 140 MeV for both models. From the analyses based 
on the wounded nucleon model, we choose the impact parameter b = 2.4 fm for 0-6 % 
central collisions and b = 8.9 fm for 35-45 % non-central collisions. The resultant number 
of participants is 342 (94) for b = 2.4 (8.9) fm, which is consistent with estimation by the 
PHOBOS Collaboration HOI. We reasonably reproduce the data in not only central but also 
non-central collisions by using initial parameters in the previous section. After tuning initial 

3 When we calculate the two-pion correlation function in the next section, we neglect this contribution for 
simplicity. 
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FIG. 6: Pseudorapidity distribution of charged particles in Au+Au 130A GeV central and semi- 



central collisions. Data from Ref. |4C]. 



parameters for central events, we have no adjustable parameters for non-central events due 
to assuming the binary collision scaling. For non-central collisions, we change the impact 
parameter in Eq. (^J) as to the number of participants with keeping other parameters in the 
initial condition. Although the binary collisions contribute to hard components, the binary 
collision scaling seems to be reasonable to parameterize the hydrodynamic initial condition 
41] from Fig. || 



We next show in Fig. |7| the p t spectrum for negative charged hadrons. The experimental 
data (0-5% central ^|| and 30-40 % semi central events H49[0 are observed by the STAR 
Collaboration. The impact parameters used in this calculation are b = 2.4 fm for central 
and b = 8.6 fm for semi central events. We represent the thermal freeze-out temperature 
dependence of p t spectrum for the models CE and PCE. The slope of p t spectrum is almost 
independent of T th near p t ~ 1 GeV/c, which is a little peculiar behavior in the usual sense. 
However, it is interpreted by the result in Fig. |5|. When one reduces the thermal freeze-out 
temperature by hand, the average radial flow enhances as its response. The magnitude of the 
response is governed by the EOS. The resultant p t slope is a competition between these two 
effects: The reduction of temperature makes the slope steeper in the case of vanishing flow, 
while the thermal distribution is Lorentz-boosted by radial flow and the p t slope becomes 
flatter. The effect of generated radial flow on the p t slope usually overcomes that of the 
reduction of T th in the model CE, so the pt slope becomes flatter as decreasing T th . On 
the other hand, the radial flow is slightly suppressed in the model PCE as shown in Fig. ||. 
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FIG. 7: Transverse momentum spectra of negative charged particles in Au+Au 130^4 GeV central 
and semi-central collisions for (a) the model CE and (b) the model PCE. The dashed, dotted, and 
solid lines correspond to T th = 100, 120, and 140 MeV. Data from Refs. EM and @. 



Hence the reduction of T th is just compensated by its response to radial flow for the model 
PCE. This is the reason why the pt slope is almost independent of T th in Fig. |7|. For the 
model CE, we reproduce the slope by choosing T th = 140 MeV. For the model PCE, in 
any T th within a plausible range, we reproduce the experimental data below 1 GeV/c for 
central collisions. This indicates that there exists the onset of hard processes around pt ~ 1 
GeV/c. It should be noted that a bend of the spectrum in low p t region is simply due to 
the Jacobian between the pseudorapidity 77 and the rapidity Y . 

We next show the transverse momentum dependence of the second Fourier coefficient for 
azimuthal distribution for the models CE and PCE. Figure M represents v^ijpt) for charged 



particles in minimum bias collisions. Experimental data is also observed by STAR [50 
Hydrodynamic analysis of the data has already been done by Kolb et al. [ pl| . They found 
the hydrodynamic result excellently coincides with the data below p t ~ 1.5 GeV/c. Our 
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FIG. 8: V2(pt) for charged hadrons in Au+Au 130A GeV collisions. The dashed, dotted, and solid 
lines correspond to T th = 100, 120, and 140 MeV. Data from Ref. |0|. 



results for the model CE are consistent with their results. Now we find the space-time 
evolution in our model is different from their result where chemical equilibrium is always 
assumed, we must check whether the hydrodynamic description is really good at RHIC even 
when we include the effects of early chemical freeze-out. The numerical results are calculated 
from the following equation: 



V2{Pt 



J2b I drjdcf) cos(20)6 



(IN 



p t dp t dr)d<f> 



(pt,r),<f);b) 



(11) 



E b fd V bJ^( Pt ,r,;b) 
Here the summation with respect to the impact parameter b is taken over every 2 fm up to 
14 fm in this analysis. The integral region of r] is from —1.3 to 1.3, which corresponds to the 
analysis by STAR |5(|. The value of v 2 (pt) depends on the thermal freeze-out temperature 
T th in contrast with the pt spectrum. We also reproduce the experimental data below 1 
GeV/c by choosing T th = 140 MeV and slightly overestimate the data above 1 GeV/c. 
Similar to the p t spectrum, this result also indicates that the hard contribution, which 
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FIG. 9: v%(ji) for charged hadrons in Au+Au 12>QA GeV collisions. The solid and dashed lines 
correspond to the models PCE and CE. Data from Refs. [50] and [53|. 

reduces v 2 calculated from hydrodynamic source H |5"2fl , becomes important above 1 GeV/c. 

In Fig. [|, we compare the pseudorapidity dependence of elliptic flow between the models 
CE and PCE. Similar to the Vi{j>t)-> ^lij]) which is to be compared with minimum bias data 



is 



v 2 {rj) 



Y^b I Ptdptdcj) cos(20)6 



(IN 



ptdptdrjdif) 



(Pt, V, <P\ b) 



(12) 



We choose T th = 140 MeV and integrate with respect to p t from to 2.0 GeV/c. Data 



plots are observed by PHOBOS |53| . The rectangular area corresponds to the statement by 
STAR H, v 2 {r]) = 4.5 ± 0.5% for 0.1 < p t < 2.0 GeV/c and | r) |< 1.3. As well as the 
case of radial flow, the elliptic flow is also reduced by taking account of the chemical freeze- 
out. We reproduce the PHOBOS data only near mid(pseudo)rapidity and overestimate in 
forward and backward rapidity. On the other hand, a microscopic transport model (JAM) 
reproduces the data only in forward and backward rapidity regions | 5"4]| . This indicates 
that the full thermalization is achieved only near midrapidity, although there are some open 
problems in hydrodynamics such as the treatment of freeze-out through the Cooper-Frye 
formula, more sophisticated initialization, and the absorption by spectators |J. 

From Figs. || and || the hydrodynamic description with early chemical freeze-out seems 
to be valid for, at least, < pt < 1 GeV/c and — 1 < r] < 1 in Au+Au collisions at 130A 
GeV. 
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FIG. 10: The transverse momentum spectra for negative pions, negative kaons, and antiprotons 
for (a) the model CE and (b) the model PCE. To see these results clearly, the yield of kaons 
(antiprotons) is scaled by 10 _1 (10 -2 ). The dashed, dotted, solid lines represent results for T th = 
100, 120, and 140 MeV, respectively. Data from Ref. @. 

B. Spectra and Flow for Identified Hadrons 



In this subsection, let us see the difference between the models CE and PCE by comparing 
Pt spectra and v^ipt) f° r identified hadrons which are supposed to be sensitive to the early 
chemical freeze out. 

The pt spectra for identified hadrons in (0-5%) central collisions observed by the PHENIX 
Collaboration [^BJ are compared with our results in Fig.|TU|. The impact parameter which 
we choose for central collisions is also b = 2.4 fm. For the model CE, the slopes of pions, 
kaons, and antiprotons become steeper as increasing T th , which is similar to the case of 
charged particles. The number of antiproton becomes very small at T th = 100 MeV for the 
model CE due to chemical equilibrium. On the other hand, the numbers of each hadron 
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FIG. 11: V2(pt) for pions, kaons, and protons in Au + Au 130^1 GeV collisions. Left (right) 
column represents the results for the model CE (PCE). Data from Ref. [|56|| . 

in the model PCE are independent of T th and reasonably agree with experimental data. 
The number of anti-proton in the model PCE might be slightly improved by taking into 
account the baryonic chemical potential which is neglected in the present analysis. It should 
be noted that pion spectra in large p t (>1.5 GeV/c) region can be reproduced by including 
contribution from non-thermalized hard partons (jets) with energy loss fl55| . 

We next show in Fig. [11] the elliptic flow for identified hadrons and its T th dependence. 
The STAR data |56] are compared with our results with or without chemical equilibrium. 
For the model CE, ^(pt) of pions is almost independent of T th , while V2{pt) of kaons and 
protons are increasing with T th . On the other hand, t>2 of pions grows with decreasing T th 
for the model PCE. Whether V2(pt) increases with T th depends not only on the particle mass 
but also on the flow velocities and its anisotropy in the transverse plane f57]. The elliptic 
flow seems to be more sensitive to T th than p t spectra when we consider the early chemical 
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FIG. 12: v^ipt) for charged pions. The solid, dotted, and dashed lines correspond to total pions, 
pions directly emitted from freeze-out hypersurface, and pions from resonance decays. Data from 



Ref. p|. 



freeze-out. For the transverse momentum spectra of identified hadron shown in Fig. [7], we 
roughly reproduce the slope in low pt region with T th = 140 MeV. On the other hand, we 
cannot reproduce V2{pt) of identified hadrons by a common thermal freeze-out temperature: 
each hadron seems to favor different T th . This indicates the hadronic afterburner in the late 
stage of the expansion may be important |[58|| . From hydrodynamic point of view, detailed 
analyses with various EOS's and initial conditions are needed in understanding thermal 
freeze-out properties. 

Figure [l2| separately represents the contribution of pions directly emitted from freeze-out 
hypersurface and the contribution of pions from resonance decays. The STAR data shows 



the contribution from identified pions [56]. From this figure, the concavity of Vzipt) m l° w Pt 
region (p t < 0.3 GeV/c) results from the resonance decay after thermal freeze-out. The exact 
treatment of the decay kinematics of resonances dilutes v 2 for direct pions especially in low 



Pt region [32|. The fraction of the contribution from resonance decays plays a very important 
role in understanding V2(pt) for pions in low pt region, so the early chemical freeze-out must 
be included when we discuss proper phenomena to the low transverse momentum. 
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VI. PION INTERFEROMETRY 



From the single particle spectra and the azimuthal asymmetry, we obtain the information 
about the distribution in the momentum space at freeze-out. On the other hand, we can 
obtain the information about the particle distribution in the coordinate space through two- 
particle interferometry [59|. We see in Sec. |I^ the space-time evolutions are considerably 



different between the models with and without chemical equilibrium. To see this more 
clearly, we discuss the two-pion correlation function in this section. 

In hydrodynamics, the two-particle correlation function for directly emitted bosons from 
freeze-out hypersurface E can be calculated from pUl pT 

P(Pl,P2) 



C 2 (Pl,P2) 



P(Pl)P(P2) 



1 + 



E dN E dN 



(13) 



Here P{px,P2) is the two particle coincidence cross section and P(pi) is the same Cooper- 
Frye formula represented in Eq. flTOl). We consider only directly emitted pions for simplicity. 

= (pi + P2V2 and g M = Pi — p^ are, respectively, the average and relative four momen- 
tum of pair, /be = (e 1 — is the Bose-Einstein distribution function. The informations 
about hydrodynamic simulations enter through the freeze-out hypersurface E and the four 
velocity in this equation. The average pair momentum K is decomposed into the trans- 
verse momentum Kt, the longitudinal momentum K z , and the azimuthal angle Qk- 4 The 
relative pair momentum g M is also decomposed into the standard coordinate, g out (parallel 
to Kt), <2Wg (along the beam direction), and g S id e (perpendicular to the others). Since the 
experimental acceptances are limited to midrapidity, | Y \< 0.5 (STAR) |62| or | r\ |< 0.35 
(PHENIX) f63|], we can put K z = 0. Moreover, we average the two-particle function C2 over 



the azimuthal angle Thus we obtain the following equation which can be compared 
with the experimental data: 

f$ K P( Pl ,p 2 ) 



CiiK^i <?side; Qout; ?long) 



f$ K P( Pl )P(p 2 ) 



(14) 

K z =0 



4 Even for central events, we have no longer cylindrical symmetry around the collision axis due to small 
(but finite) value of the impact parameter, so C 2 depends on the azimuthal angle <frj< which is measured 
from the reaction plane. 



22 



Our definition of the HBT radii is similar to the one in Refs. [61, B3, |65[. Assuming that 



C-i{Kt, q S ide, 9outj ftong) for each Kt is fitted by the Gaussian form 

C 2 = 1 + A exp(— R aide q side - R 2 out ql ut - -Rfongftong) 

with a chaoticity A = 1, the Kt dependence of HBT radius for the side direction is 
-Rside(-Kr) = V9side(-^r)j where C 2 (K T , q* ide , 0, 0) = 1 + e -1 , and analogous definitions for 
the Kt dependence of R OVLt and .Riong- 

We evaluate two-pion correlation functions for negative pions directly emitted from freeze- 



out hypersurface and obtain the K T dependence of -R S id e , -Rout, and R\ ong . Figure |1^ shows 
the HBT radii and the ratio idt/i^ido for the models CE and PCE with T th = 140 MeV. 
Here the impact parameter which we choose in this analysis is b = 2.4 fm. The difference 
between the two models is very small for -R s ;d c - On the other hand, R OVLt and -Riong for the 
model PCE are significantly smaller than the ones for the model CE, which reflects the 
space-time evolution of freeze-out hypersurface depicted in Fig. £2[ We compare our results 
with experimental data observed by STAR (12% most central events) [Q and PHENIX 
(30% most central events) p5| . We reproduce the K T dependence of R\ ong by employing 
the model PCE, while -R S ide and i? Q ut do not show good coincidences with the experimental 
data. 

It has been suggested that the ratio R ut/R S ide reflects the prolongation of the lifetime 
due to the phase transition between the QGP phase and the hadron phase [ |BT| . Various 
models predict this ratio has a value significantly larger than unity in some Kt region 
61], |64], |65], |H, although it is around unity in 0.2 < K T < 0.6 GeV/c in Au + Au 



y/sNN = 130A GeV collisions according to the recent measurement at RHIC |62|, |6~3T| . This 



discrepancy is often called "the HBT puzzle" [|6"6| , |67 [. The ratio for the model PCE reduces 
by about 12 % above 0.2 GeV/c due to the property of the early chemical freeze-out, but it 
clearly turns out to be larger than the experimental data. Our aim in this paper is to see 
how the early chemical freeze-out affects the HBT radii, so we leave detailed discussions on 
the HBT puzzle in the future works. 5 



5 In Ref. p^| , slightly larger i? s ;do is obtained by assuming simple flat transverse profile and vanishing 
impact parameter for initial energy density. The flat profile leads a larger value of variance of a fluid in 
the transverse plane than the present initialization based on the binary collision scaling. 
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VII. SUMMARY AND DISCUSSION 



We have investigated the effect of early chemical freeze-out on the hydro dynamic evolu- 
tion and the particle spectra by using a genuine three dimensional hydrodynamic model. We 
constructed the equation of state for hadronic matter which is in partial chemical equilib- 
rium. Pressure as a function of energy density is not affected by the chemical non-equilibrium 
property, while temperature as a function of energy density is largely reduced due to the 
large population of resonance particles. By using the EOS with a first order phase transi- 
tion, we simulated hydrodynamic evolution at the RHIC energy. We found the system cools 
down more rapidly than the conventional model and that the lifetime and the spatial size of 
the fluid and the radial and elliptic flows are reduced when we compared at isothermal hy- 
persurface. We also analyzed particle spectra and two-pion correlation functions in Au+Au 
130A GeV collisions. We chose initial parameters in the hydrodynamic simulations so as to 
reproduce the pseudorapidity distribution in central collisions observed by PHOBOS. The 
slope of p t spectrum for negative hadrons is less sensitive to the thermal freeze-out temper- 
ature, which results from the reduction of radial flow. On the other hand, the transverse 
momentum dependence of elliptic flow v 2 (pt) for charged hadrons depends on the thermal 
freeze-out temperature. The situation is completely opposite to the ordinary hydrodynamic 
results. In the conventional models, the p t slope is steeper as increasing T th and i>2(pt) 
is not so sensitive to T th . We found T th = 140 MeV and the resultant average radial flow 
< v r >= 0.38c are the values to simultaneously reproduce the p t slope and i>2(pt) for charged 
hadrons (mainly pions) below 1 GeV. 1*2(77) is also reduced by the early chemical freeze-out, 
but we failed to reproduce the data in forward and backward rapidity regions. We see the 
thermal freeze-out temperature independence of pt slope more clearly in the pt spectra for 
identified hadrons. We reasonably reproduce the p t slope with T th = 100-140 MeV for the 
model PCE, while i>2(pt) for identified hadrons seems to favor different thermal freeze-out 
temperatures which are dependent on hadronic species. In order to see more quantitatively 
the effect of early chemical freeze-out on the temporal or spatial size of fluids, we calculated 
the two-pion correlation functions and obtained the HBT radii. By taking account of the 
early chemical freeze-out in the EOS, -R out , -Riong, and R out /R S idc are significantly suppressed, 
while i? S ide is not changed. Nevertheless, the properties of the early chemical freeze-out are 
not enough to interpret the HBT puzzle. 
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It would be very interesting to see other observables in the case of partial chemical 
equilibrium. Although the best place to see the difference between the models CE and PCE 
must be the particle ratios, we cannot discuss this observable by using the present model 
due to an approximation of vanishing baryonic chemical potential. This will be discussed 
elsewhere. Penetrating probes such as thermal photons and dileptons may also be affected 
by the early chemical freeze-out. The emission rates of photons or dileptons increase due 
to the chemical potentials for hadrons, while the space-time volume of the hadron phase 
decreases. Therefore we should check whether the total multiplicity and spectra of photons 
and dileptons are changed in terms of hydrodynamics. 

There are similar dynamical approaches to describe the early chemical freeze-out. The 
model evolves the QGP and mixed phases as a relativistic fluid, while it switches to a 
hadronic cascade model, e.g., UrQMD [Q or RQMD [f)8[j . The advantage of our hydrody- 
namic model over these hybrid (hydrodynamics + cascade) models is to be able to obtain 
the average hydrodynamic behavior and the effect of temperature naturally. When one dis- 
cusses the spectral changes of hadrons due to the medium effects such as temperature and/or 
baryonic chemical potential [E5], one can easily estimate its effect by using hydrodynamic 
simulations. We regard the model PCE as a complementary tool to those hybrid models. 

After publishing our preliminary results and almost finishing this work, the authors 
are aware of a paper concerning the same subject [J70], in which the conclusion on the 
hydrodynamic behavior is almost the same as ours. In Ref. J7[J, the effect of finite baryonic 
chemical potential is included but the hydrodynamic simulation is performed only in the 
transverse plane by assuming the Bjorken's scaling solution. 
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APPENDIX A: MONTE CARLO CALCULATION OF THE CONTRIBUTION 
FROM RESONANCE DECAYS 



In this appendix, we show how to calculate the particle distribution from resonance decays 
within the Cooper-Frye prescription. This method can be used in hydrodynamics-motivated 
models as well as in hydrodynamics. 

Lorentz transformation for the momentum of a decay particle between the local rest 
system (starred) and the finite momentum system (non-starred) of a resonance particle R is 



P =P~PR 
We rewrite Eq. ( |AT| ) explicitly 

Pi 



E 



PPr 



m R m R (m R + E R ) 



cos ( 



where 



F(p 



Pi -pmF(pi,(j>), 
Pi 
Pi 

Ptipu 4>) cos0 - p Rt cos (ft R F(pi, 0) 

y/Ptipu 4>) + P 2 Rt F2 iPh </>) - 2pt(pi, (t>)PRt cos(0 - <p R )F(p h 

E{pu 0) Pt{Ph </>)Pm cos(0 - (j) R ) + pipm 



(Al) 



(A2) 



(A3) 



(A4) 



m R m R (m R + E R ) 

Here the independent variables which we choose for decay particles are the longitudinal 
momentum pi and the azimuthal angle <fi. Thus the transverse momentum of a decay particle 
Pt is written in terms of pi and 

Pt(pi,4>) = f, 2 1 2 , t — -tt\(E* +p l v R aR)v Rt cos( y (p- R ) 

± \j (E* + piv R a R ) 2 - (pf + m 2 )7| [1 - v 2 Rt cos 2 (0 - <f> R )] \ . 
The Jacobian of the Lorentz transformation is defined by 

dp*dcj)* = J(p h <f>;V R )dpid<f>, 



(A5) 



APi 



R) 



(A6) 
(A7) 



dpi d<t> 
d<t>* d<t>* 
dpi dcj> 

The calculation of J is straightforward, so that we do not represent it here. The normaliza- 
tion of momentum space volume for a decay particle in the resonance rest frame is 



*>* dg_ 
P * 2p* 



2t 



2tt 



1. 



(A8) 
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We always average the decay probability over the spin of resonances, so that the decay 
probability does not depend on pi and 0*. Thus the normalization in the resonance reference 
frame is 

J(piA; VR)dpid(j) = 1 

The Jacobian in Eq. (|A7[ ) has very narrow peaks when the resonance particle moves at a 
large velocity in the laboratory system |fffi| . This singularity makes it difficult to integrate 
the Jacobian numerically. So we introduce a very simple Monte Carlo calculation to evaluate 
the momentum distribution from resonance decays. All input parameters in this calculation 
are the numerical results of hydrodynamic simulation, i.e., the temperature T th , the chemical 
potential for resonance particles /zr, the three-dimensional fluid velocity v, and the element 
da^ on the freeze-out hypersurface S. In the following discussion, we show how to obtain 
the rapidity distribution of negative pions, for simplicity, only from p mesons. In this case, 
the branching ratio B p o(-)^ n - n +(o) = 1. It is straightforward to extend this scheme to the 
cases for other resonances or the transverse mass (momentum) distribution. 
Step 1 : Evaluate the number of p° and p~ which are emitted from or absorbed by the k-th 
freeze-out hypersurface element da^: 

t r _ 9R f d 3 p R | p R ■ da k 1 



k (2tt)3 J s E R exp [(p R ■ u k - p R ) /T th ] - 1 ' 1 } 

The integrand does not contain the Jacobian, so that it is simple to carry out the numerical 
integration by a standard technique. It should be noted that Nj? is different from the net 
number of emitted p mesons from the k-th. fluid element and that this is used merely for 
normalization. 

Step 2: Generate random momenta P* (1 < j < N) for p mesons which obey the 
distribution 

p*2 

(All) 



exp 



(y P* 2 + ml - fij^ /T th 



- 1 



Here we omit the Breit-Wigner function for simplicity. 

Step 3: For each A^ random momentum P*, generate random variables (G*,$*) whose 
ensemble is uniformly distributed on the unit sphere. By using these random vari- 
ables, we obtain an ensemble of p meson with momentum P* = (P*j,P*j,P*j) = 
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(P* sin 6* cos P* sin 0* sin $*, P* cos ©*), which obeys the Bose-Einstein distribution in 
the fluid rest system. 

Step 4 '■ Boost P* with respect to the fluid velocity v k 

* = ^™fa + ^)' (A12) 

where 7 fc = - v\. 

Step 5: Generate N uniform random variables on the unit sphere (9*,(j)j) and 

obtain an ensemble of negative pions with momentum p* = {Pxj,Plj,P* Z j) — 
(p* sin 9* cos <f>j, p* sin 9* sin (f)*,p* cos 9*) , where p* is given by 

P* = ^~\/ ( m R + m ^) 2 ~ m x\J ( m R ~ m n) 2 - rn? x . (A13) 
For the decay process p — > 7r7r, = m^. 

S'tep Boost p* with respect to the resonance momentum Pj 



Pj = P*j + p j 
Step 7: If Pjda^k is positive, 



E 3 + P, i\ 



m R m R (m R + E R ) 



(A14) 



K - K + (A15) 



If P^da^k is negative, 



«r - jv? + ^^p^. (aw) 

Here, N£ (N^) is to be proportional to the number of p mesons which are emitted from 
(absorbed by) the k-th fluid element. 

Step 8: If the rapidity of a negative pion Yj which is evaluated from pj enters in a rapidity 



window Y — 4p < < F + 4p and P^da^k is positive, 



AiV+(F) - AiV+(F) + (A17) 



If Yj also enters the above rapidity window but Pjda^k is negative, 



AN~(Y) -> AiV-(F) + |Pj ^ fcl . (A18) 
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Step 9: Repeat steps 7 and 8 for all N random variables. 

Step 10: Obtain the rapidity distribution of decay particles from the k-th. fluid element 



dY { ' N+ + NZ AY 



(A19) 



k ~ k I 

It should be noted that the minus sign in the bracket of Eq. ( |A19|) means the net number 
of emitted particles from the k-th fluid element, which is consistent with the Cooper- Frye 
prescription. For the normalization in Eq. (|A19| ), we use the gross number N£ + N% since 
this number is positive definite. 

Step 11: Repeat the above steps from 1 to 10 for all fluid elements obtained in a numer- 
ical simulation of the hydrodynamic model. Summing over the contribution from all fluid 
elements on the freeze-out hypersurface S, we obtain the rapidity distribution of negative 
pions which are from p decays: 

' B SF L ( Y ) = j: ! W(n (A20) 

k 
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